RECONSTRUCTION FROM RADON PROJECTIONS AND 
ORTHOGONAL EXPANSION ON A BALL 
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Abstract. The relation between Radon transform and orthogonal expansions 
of a function on the unit ball in R'' is exploited. A compact formula for the 
partial sums of the expansion is given in terms of the Radon transform, which 
leads to algorithms for image reconstruction from Radon data. The relation 
between orthogonal expansion and the singular value decomposition of the 
Radon transform is also exploited. 



1. Introduction 

Reconstruction of an image from its Radon projections is the central theme in 
x-ray tomography and has spectacular applications in medical imaging. Mathemat- 
ically the problem is to find a good approximation to a function based on a finite 
collection of its Radon projections (see, for example, [9, 10, 19]). 

The main topic of this paper is the connection between the Radon transform 
and the orthogonal expansion of the function on a unit ball. This connection 
was initiated in the classical paper [4] with an inversion formula of the Radon 
transform based on spherical harmonic expansions. The relation between the Radon 
transform of a function, supported on the unit ball, and its orthogonal expansion 
was studied or used in [5, 6, 11, 12, 15, 17], among others (see [19] for further 
references). The papers [5, 6, 12] studied also the singular value decomposition 
(SVD) of the Radon transform using an orthogonal basis. Since then SVD has 
become an important tool for studying the stability of the inversion problem, the 
resolution of the reconstruction, and the incomplete data problem; see, for example, 
[3, 6, 13, 14, 19]. The truncated SVD also provides an algorithm for reconstruction 
of images. Because of the complicated formulas involved in the orthogonal or SVD 
expansions (see, for example, [5, 12, 19]), the algorithms did not seem to be used 
in practical applications. 

Recently a new reconstruction algorithm was proposed in [27] and further stud- 
ied in [28, 29]. The new algorithm is called OPED, as it is based on orthogonal 
polynomial expansion on the unit disk. The algorithm reproduces polynomials of 
high degrees and allows a fast implementation ([28]). The numerical tests shows 
that the algorithm is fast, stable, and produces high quality images ([28, 29]. The 
key ingredient for deriving the algorithm is the following formula for the partial 
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sum S2mf of the orthogonal expansion of / on the unit disk, 

(1.1) S2mf{x,y) = - — — tX] / T^'I'u f{t)^'2m{t,x COS (j)^ +y sin (l)^)dt, 

where </)j, = 2m+i ^^'^T^ofit) is the Radon projection on the hne x cos sin = t 
(see Section 3). It turns out that there is a natural extension of this formula to 
the unit ball of higher dimension, which shows that the orthogonal polynomial 
expansion of / can be expressed in terms of the Radon transforms and allows us to 
extend the OPED algorithm in the unit ball of W^. Furthermore, there is a close 
relation between SVD and the extension of the formula (1.1). In fact, they can 
be brought together by the use of a compact formula of the reproducing kernel of 
orthogonal polynomials in [24] . The orthogonal expansion on the unit ball has been 
studied recently in [24, 26], which can be used, in particular, to derive the uniform 
convergence of the algorithms. 

The purpose of this paper is two folds. Firstly we will clarify the relation be- 
tween orthogonal expansion on the ball and the Radon projections and derive the 
extension of the OPED algorithm in higher dimensions. Secondly, we will explain 
the connection between SVD of the Radon transform and orthogonal expansions. 
In particular, wc shall show that using truncated SVD to reconstruct the image is 
the same as using OPED algorithm. 

The paper is organized as follows. The following section contains a succinct 
account of the basic results on orthogonal polynomials on the unit ball. The or- 
thogonal expansions in terms of the Radon projections is developed in Section 3. 
The extension of the OPED algorithms and a convergence result are given in Section 
4. Finally, the SVD of the Radon transform is discussed in Section 5. 



2. Preliminaries on orthogonal polynomials 

Let B"^ := {x : \\x\\ < 1} and 5'*"^ := {a; : ||a;|| = 1} be the unit ball and the 
unit sphere of R'', respectively. We denote the surface area of S'^~^ by aj, and the 
volume of B"^ by bd- Then 

= , and Od = — 



T{d/2) " d r((rf + 2)/2)' 

Inner product on the ball. For later discussion let us introduce a weight function 
on the unit ball, 

W^,{x) = (1 - \\xff-^/^, X e B'^. 
The inner product on the unit ball is denoted by 

{f,9)L^Bd) = at, f{x)g{x)Wf,{x)dx 

where is the normalization constant of W^, that is, = 1/ J^a W^{x)dx. For 
/Lt = 1/2, it is the unit weight (Lebesgue measure) and is equal to b^^. We will 

mainly work with the Lebesgue measure, so the inner product {f,g)L^[B'i) should 
be regarded as with /i = 1/2 unless specified otherwise. 

Polynomial spaces. Let IIJ^ denote the space of polynomials of degree n in d 
variables. We say that P € 11^ is an orthogonal polynomial on B'^ if (P, Q)l'^(^B'1) = 



RADON PROJECTIONS AND ORTHOGONAL EXPANSIONS 



3 



for all Q S ^n-i- Let denote the space of orthogonal polynomials. It is well- 
known that 

dimnf = and dimV^ = ' 



n J \ n 

Several explicit orthonormal bases of are known (see, for example, [7]). We will 
need one given in terms of the Jacobi polynomials and spherical harmonics. 

Jacobi polynomials. The fc-th Jacobi polynomial is dentoed by pI"'^^ and they 
satisfy the orthogonal relation ([23]) 

(2.1) c„,^ /' Pt'-^^it)P^"^''^(t)w^,p(t)dt 



1 



(a + l)fe(/3 + l)k{a + /3 + fc + 1)^^ ^ _^ h^^'^^Sk i 



k,, 



fc!(a + /3 + 2)fe(a + /3 + 2/c + 1) 
where Wa,f3{t) = (1 — + '^a,/3 is the normalization constant of Wa,f3, 

and the notation {a)k ■= a{a + 1) • • • (a + fc — 1) denotes the shifted factorial 
(Pochhammer symbol). From (2.1) the orthonormal Jacobi polynomials are given 

hy p^^'^\f) := [/i(r'«]-i/2p^/5)(i). 

Gegenbauer polynomials and Chebyshev polynomials. When a = P = X — 1/2, 
the Jacobi polynomials become the Gegenbauer polynomials, usually denoted by 
Cj^ and normalized by 

(2.2) c. |'^C,^t)Q\t)(l-t^)^-i/^dt=^^|^<5.,:=4^)5, 

where cx = r(l/2)r(A + l/2)/r(A + 1). When A = 1 and A = 0, Cj^{t) becomes 
the Chebyshev polynomial of the second kind, Uk{t), and the first kind, Tk{t), 
respectively, and 

sinfA; -I- 1)9 

(2.3) Uk{t) = — ^^.^^ ' and Tk{t)=cosk0, where t = cose. 

Spherical harmonics. These are defined as the restriction of the homogeneous 
harmonic polynomials on the sphere. Let H!^ denote the space of spherical har- 
monics of degree n in d variables. It is known that 

^n + d-l\ /n + d-3'^ 
n J \ n 

Let {Yk^n : 1 < < dimW^} denote an orthonormal basis of 'Hf^. Then 
<^d^ [ Yk,n{OYi,n{OMO = l<k,l<dimnt 



dimH^ 



d 



We emphasis that Yk^n{x) arc in fact homogeneous polynomials in II! 

An orthonormal basis for V^. We give the basis for inner product defined in 
terms of W^(a;). Setting fj. = 1/2 gives the basis for the Lebesgue measure. Let 
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Yj,m be an orthonormal basis for TiJ^^. Define 

(2.4) f}:^,{x) = [hnA-'Pk~'''"~"'^'^\nAf - m,n-2k{x), 

where 

_ r(/x + ^)r(n-2fc + f) 
[tinM -r(|)r(n-2fc + M + ^)" 

Then the set {f^j : I < j < dim W^_2;,, < 2k < n} is an orthonormal basis for 
V^; that is, e and {fljJ^,j,)L-iB^) = 5k,k'hj' (see [7, p. 39]). 

Reproducing kernel ofV^- The reproducing kernel Pn{-, •) of satisfies 



(2.5) f{y)P„ix,y)W^iy)dy ^ fix), V/ G 

Let {P;^ : 1 < A: < dimV^} denote any orthonormal basis of V^. Then 

iV„ 

Pn{x, y)=J2 Pk{^)Pk(y), = dim V^. 

fe=i 

The definition of P„(-, •), however, is independent of the particular choice of bases. 
In particular, we can take the orthonormal basis in (2.6) and get 

0<2k<n j=l 

The reproducing kernel satisfies a compact formula that will play a fundamental 
role in our study; it is given by ([24]) 

Pn{x,y) = ^c^-i £^ C,>:i{x,y) + ^1 - Wxr^/T^^ s){l - s^-'ds 

where A = /x + (., •) is the Euclidean inner product in M*^, and cx is defined in 
(2.2). In particular, it implies that 

(2.7) p^ix,0 = '^C^{{x,0), ^&S''-\ xeB''. 

Orthogonal expansions on B'^ . If {PJ} : 1 < fc < iV„}, 7V„ = dimV^, is an 
orthonormal basis of V^, then the standard Hilbert space theory states that there 
is an orthogonal expansion 

OO Nn 

k=0 k=l 

Let proj^ : L'^{B'^) i-^ denote the projection operator. Using the reproducing 
kernel, the orthogonal expansion can be stated as 

OO ^ 

(2.8) / = EP™jA/' "^^^^^ Projfc/ = a^/ f{y)Pn{x,y)Wi_,{y)dy, 
which is independent of the particular choices of the bases of V^. 
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3. Radon Transform and Orthogonal Polynomial Expansion 

Let / e -L^ be a real valued function. For ^ e S"*"^ and t £ R, the Radon 
transform of / is defined as 

nm,t)--= [ fix)dx= [ f{t^ + y)dy, 

where the integral is over a hyperplane of {d — l)-dimension perpendicular to ^ and 
with minimum distance t to the origin. More general definition on other spaces 

or manifolds can be found in [8] . For properties of Radon transforms wc refer to 
[8, 19]. We assume that / has compact support in B'^, so that the integral above 
should be understood as over B'^ f\{x : {£,,x) = t}. In particular, for ^ e S'^~^, let 
denote an orthogonal matrix whose first row is ^ and let B'^{r) denote the ball 
of radius r in M'^; then a change of variables x i-^ {t, y)Q^ shows that 

(3.1) nfiu)= [ f{it,y)Qi)dy 

= {i-t^)^ f ,f{{t,VT~^y)Qddy. 

Since {{t,y)Q^,S,) = t, an immediate consequence of (3.1) is the following identity, 

(3.2) / f{x)g{{x,0)dx= [\f{^,t)g{t)dt, ^ G S''-\ 

JB'i J-1 

whenever both integrals make sense. The definition of TZf also implies that 

(3.3) nf{-^,-t) = nf{tt), ^gS''-\ tGR. 

For fixed ^ and t, wc also call Tlf{^, t) a Radon projection. The essential problem 
for x-ray imaging is to find a good approximation to the function / based on a 
given data set of its Rdaon projections. 

Wc now derive the orthogonal expansion of / in terms of Radon projections. 
The following proposition plays a key role. 

Proposition 3.1. For x,y £ B'^ , the reproducing kernel Pn(-, •) satisfies 

(3.4) P„(x,y) = ^^CT-' I Ct'\{x,^))Ct/^{{yS)d^{0. 

Proof. Prom the explicit formula of fj^ - at (2.4) with = 1/2, we deduce that 

(3.5) fl^iO = HnY^,n-2m, ^eS''-\ 

where, using the fact that pf'^\t) = [h^'^Y^/^Ph'^^Ht), -Pf '^'(1) = 1, and the 
formula of 4"''^' in (2-1), we have 

/q rj r. 1-1 (0,n-2fe+^),,^ n + d/2 

(3.6) Hn = [hn,k\ Pk (1) = d , 

independent of k. Consequently, integrating over S'^~^ we get 

Multiplying the above equation by fkj{x) and f^ijiiy) and summing over all 
k, k', the stated result follows from (2.6) and (2.7). □ 
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Theorem 3.2. For n>0, 

In particular, for f e L'^{B'^), 

(3-7) / = E ^ir<^a' I C nf{U)C'J'{t)dtC'J\{;0)du^{0- 

^0 '^/^ Js^-^ J-1 

Proof. By the formula (2.5) with /x = 1/2 and the formula (3.4) of P„(-, •) we have 

Jb<^ 

= ^^^"^^ I / f{y)ct'\{mmyct'\{x,0)d^{^). 

The identity (3.2) shows that the inner integral is 

(3.8) / f{y)Ctl\{y,i))dy= f nf{i,t)Ct/\t)dt, 

SO that the stated formula follows. □ 

The formula (3.7) as stated here has already appeared in [20] in a study of the 

approximation by ridge functions. See also [1] for the case of d = 2. Although 
spherical harmonics expansions for d = 2 was used in the classical work of [4], 
its compact form in (3.7) is quite recent and not used for reconstructing images 
from Radon data until recently ([27]). It should also be noted that for li > 2, the 
Gegenbauer polynomials and spherical harmonics were used for constructing Radon 
transforms already in [12]. 

Let us mention that there does not seem to be an analogOTis formiila for the more 
general case of orthogonal expansion with respect to W^. In fact, in the general 
case, the formula (2.4) gives 

where 

0, rr (m + l/2)fc(M + ^)n-k{n + M + ^) 



fc!(|)„_,(M + ^) 

which depends on both n and k (comparing with (3.6)), so that Proposition 3.1 

with Cn^ replaced by Cn^ ^ does not hold. 

Let Snf denote the partial sum operator of the orthogonal expansion (2.8), 

n 

(3.10) 5„/(a;) = Eprojfc/(a;). 

Evidently, the expansion (2.8) holds in the sense that Snf ^ / in LF'{B'^) norm. 
Corollary 3.3. Let Sn be the partial sum operator defined in (3.10). Then 

(3.11) Snfix) = / nfi^, t)^n{t, {x, 0)dtduj{0. 

Jsd--^ J-1 
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where 

k + d/2 d/2, ■,^d/2 



(3.12) Mt,u) ■.= j:^^ct^\t)Ct^\u). 

A cubature formula on of degree M is a discrete sum such that 

N 

(3.13) a^' / m)MO = E ^-/(^-)' / e ^MiS"-'), 

where IiM{S''-~^) is the space of spherical polynomials, that is, the space of Hj^ 
restricted on 5'^"^. If all are positive, the cubature is called positive. We call 
a polynomial P € 11^ even if it satisfies P{x) = P{—x) for all x € M''. The 
cubature formula (3.13) is called symmetric, if it is exact for all even polynomials 
innM(5<^-^). 

Proposition 3.4. Suppose (3.13) is a symmetric cubature formula of degree 2n. 
Then 

(3.14) Snf{x) = A.&d ' / nf{^^,t)^n{t, {X, ^.))dt. 

i/=l "'-I 

Proof. The equation (3.8) shows that Px{0 '■= j\T^fiS,,t)^ni^,t',x)dt is a poly- 
nomial of degree at most 2n in ^. Furthermore, using the fact that TZfi—^, —t) = 
'R-f{^,t), it is easy to see that Px is even, so that the cubature formula on S"^"^ is 
exact when applied to Px(0- ^ 

We consider some special cases of lower dimensions below. 

The case d=2. For ^ G 5^ we write ^ = (cos 61, sin (^) and we shall write TZgf{t), 
9 G [0,27r], instead ofTZf{^,t). Since 62 = tt and the following cubature formula 

i- J^^ fiOd^iO = E mi = (COS sin ^) 

is symmetric and of degree 2n, we conclude that 

(3.15) Snf{x) = —-Y ng J {t)^n{t,Xi cos e^ + X2 sin e^)dt 

where 6*,, = ^ and is (3.12) for d = 2, 

n 

^n{t,u) = Y.^k + l)Uk{t)Uk{u). 

k=0 

This formula can be found implicitly in [11] (see (5.9), (4.3) and (3.7) there). In 
the case of n = 2m, we can use the elementary relations 

cos ^ = - COS itelz: , sin = - sin 

and the fact that TZ{6 + tt, —t) = TZ{6,t) to rewrite (3.15) as 

(3.16) S2mf{x) = - — ttE / T^'P^f {t)^2m{t, xi COS 4)1, + xi sin (j)u)dt 
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where = This is the formula (1.1) proved in [27] from which the OPED 

algorithms are derived. □ 

The case d=3. For ^ e 5^ we use the spherical coordinate 

^ = (sin (f) sin 9, sin (f) cos 9, cos (f)), < 4> < t^,0 < < 

Several explicit cubaturc formulas on the sphere arc known, sec. for example, [18, 
21]. Let tk = cos6k, k = 0,1, ... ,n, denote the zeros of the Legcndre polynomial 
of degree n + 1 and be the corresponding weights of the Legendre- Gaussian 
quadrature formula. Let 

£,k,u = {sm-^sm9k,cos-^sm9k,cos0k), Q<k,v<n 

Then the product type cubature formula 

■in ^ n n 

fc=0 i/=0 

is symmetric and of degree 2n. Consequently, we have 

(3.17) 5'„/(x) = — -^Afc^ / nf{^k,.,mn{t,{^k,.,x))dt, 

where is the function (3.12) for d = 3. □ 

The formula of Snf in terms of Radon projections allows us to give an approxi- 
mation to / based on finite Radon projections. The convergence of Snf to / holds 
in norm but does not hold in the uniform norm in general. In fact, it is known 
that [25] 

(3.18) ||5„||co = 0(n^), d>2, 

where || • ||oo is the operator norm of Sn in C{B'^), and An = 0{Bn) means ci^„ < 
Bn < CiAn for two constants Ci and C2 independent of n. There is, however, a 
simple construction that gives a better convergence result. 

Let r? be a C^+2(M) function such that 7?(t) > 0, ri{t) = 1 for < t < 1 and 77 
has compact support on [0,2]. Define 



(3.19) Slf{x):=Y,ri[-]wo]kf{x). 

fe=o 

The operator satisfies the following properties [26]: 

Proposition 3.5. Let f e LP{B'^), 1 < p < 00 or f G C{B'^) for p = 00. Then 

(1) 5;?/ = / z//en„.- 

(2) forneN, |lr,„/||p < c||/||j, 

(3) for neN, \\f - r]n\\p < cEn{f)p := inf^en^ ||/ - p\\p. 

As S^lf is a polynomial of degree 2n, the last property shows that, up to a 
constant multiple, it is close to the polynomial of the best approximation to /. 
Since proj„ / can be written in terms of Radon projections, so can S^f. 
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4. OPED ALGORITHMS FOR RECONSTRUCTION OF IMAGES 

The essential problem in computerized tomography is to find a good approxi- 
mation to the function / based on a set of discrete Radon data. The expression 
(3.14) allows us to derive such an approximation by a simple quadrature formula 
on [—1, 1]. Because of (3.1), we choose the quadrature formula to be of the form 

-1 j=0 

where 0^/2 is defined as in (2.2), and assume that it is exact for polynomials of degree 
M. In particular, we can choose the Gaussian quadrature, for which tj = tj^n, 
< i < JT., are zeros of the Gegenbauer polynomial C'^/-^i{t) and Wj are all positive 
and given by explicit formula (see [23]). The Gaussian quadrature formula is exact 
for polynomials of degree up to 2n + 1. 

Proposition 4.1. Let (3.13) be a positive symmetric cubature formula of degree 
2n and (4.1) he the Gaussian quadrature formula. Define 

N n 

(4.2) Anf{x) = 6,l^A^^w;,7^/(e.,^,)$„(^„(x,U)■ 

i/=l j=0 

Then Anf preserves polynomials of degree n, that is, Anf = f whenever f G V^. 

Proof Wc start from (3.14). If / is a polynomial of degree at most n then, by (3.1), 
(1 — t^) ~'R-f{£,vi t) is a polynomial of degree n. As $n(i, {x, ^u)) is a polynomial 
of degree n in t and the Gaussian quadrature formula is of degree 2n + 1, the fact 
that Anf = f follows. □ 

The functions Anf are obtained from the orthogonal partial sums Snf of f by 
applying the Gaussian quadrature formula. They provide a sequence of approxi- 
mation to / based on the sot of discrete Radon data 

{7^/(e.,^J): l<i^<N,0<j<n}. 

In other word, An provides an algorithm for reconstruction of images from the 
Radon data. We will show that Anf converges to / uniformly if / is smooth 
enough. First, however, we consider some special cases. 

The case d=2. In this case we can start from the formula of S2m at (3.16). The 
Gaussian quadrature formula is 

-. „1 , 2m 

- / /(t) Vl - f^dt = y sin^ V,/(cos ej), Oj = rr^^, 

■kJ_i^' 2m+l^ ^'^^ ^' ^ 2m +1 

which leads to the OPED algorithm of type II, 

2m 2m 

(4.3) Aimf{x) = J2 J2t^<^ J {cos ej)Tj,4x), 

i/=Oi=l 

where 

^ 2m 

Tj^x) = — — ——2 V](A: + l)sin((A; + l)0j)Uk{xi cos^j. + X2 sinc/)^). 
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The OPED of type II is closely related to an algorithm in [2] , where the connection 
to orthogonal polynomial expansion was not considered. The formation of the lines 
on which the Radon projections take place is often refereed to as scanning geometry, 
as it determines how the object being examined is scanned by the x-rays. We can 
use the Gaussian quadrature formula for the Chebyshev polynomials of the first 
kind, 



TT 



k—O 



to discretize the integral in (3.16) by applying it to the integrant multiplied by 
1 — leading to the OPED algorithm of type I with a different scanning geometry, 
which has the same formula as (4.3) except that 9j need to be replaced by ijjj and 
the summation on j starts from j = 0. We refer to [29] for the discussions of these 
two scanning geometries and their implementation in practical problems. 

Both types of these two OPED algorithms work well in our numerical testing ([28, 
29]). It should be mentioned that the explicit formula of Unit) in (2.3) permits a 
fast implementation of the OPED algorithm, which uses fast Fourier sine transform 
and an interpolation step ([28]). □ 

The case d = 3. In this case we can start from the formula of 5„/ at (3.17). We 
apply the Gaussian quadrature formula 



o ^1 n 

- / fm-e)dt = Y,^^f{tj), 



i=o 

where tj, < j < n, are zeors of C^/^(<). We can also apply the Gaussian quadra- 
ture formula for the Lebesgue measure. This leads to a three dimensional OPED 
algorithm, 

n n n 

(4.4) Anfix) = — y ^Afe^^w;j-7e/(a,^,ij)$„(tj-,(a,>.,a;)). 

fc=0 1^=0 j=0 

The Radon data used in (4.4) are integrals over planes {x,^k,iy) = tj. Such data 
can be approximated by integrals over lines. □ 

For d = 3, one can uses multiple 2D slices to reconstruct image on a cylindrical 
domain, as proposed in [27]. An interesting question is to see which of these two 
algorithms are more suitable for the 3D reconstruction. 

Next we consider the convergence of Anf in the uniform norm on B'^. 

Theorem 4.2. The uniform norm of the operator An is given by 

N n 

(4.5) ll^lloo = sup A„(x), A„(x) = VA^Vw;j(l-tj2)T^|i>„(tj,(a;,C^))|. 

Furthermore, there is a constant c independent of n, such that 

(4.6) II^IU < cn^''. 

In particular, if f is smooth enough then A2nf converges to f uniformly on B'^. 
Proof. To estimate the norm of An, we first observe that 



[i-t'r^nfi^^^t) 



< bd- 



RADON PROJECTIONS AND ORTHOGONAL EXPANSIONS 



11 



from which it follows that 

N n 

\\Anf Woo < ll/lloo E E "^^(1 - ^i)"^ ' 
j/=l j=0 

since bd-ib'^^ — Cd/2- Taking the maximum over B'^ shows that ||^|joo is bounded 
by the right hand side of (4.5). To prove the equal sign, we construct a function 
fs for each e > such that ||/e||oo = 1 and H^/sHco > max^g^d A„(x) — ce. A 
moment of reflection shows that the construction can be carried out easily; see [27] 
for one special case of d = 2. 

To prove (4.6) we use (3.12) and the fact that \C^{t)\ < C^(l) = = 
0{n^^^^), which implies that 

< E '-i^iC'J'il)? < c± '-i^k^'-' < on-'. 

k=0 ' k=0 ' 

Since A^^ and wj are all positive and, as the cubature and the quadrature are 

exactly for constant function, X^^j = 1 and J=o = 1, we conclude that 
ll^nll < cn^''. If / e C^"*, then the fact that A„p = piov p€U.^ and the triangle 
inequality shows that 

\\Anf - /lloo < (1 + \\A\\oo)EM)oo < cn'"'EM)oo- 

It is shown in [27] that En{f) < cn"^'' ||'D''/||, where D is a second order differential 
operator, so that the convergence of Anf for functions smooth enough follows. □ 

We should point out that the estimate (4.6) is a rough upper bound, the actual 
norm should be smaller. In fact, in the case of d = 2, the norm of A2m at (4.3) was 
estimated in [27] to be 

IIAmlloo ^ mlog(m + 1), 
which is sharp and is just slightly worse than the estimate (3.18) of the norm of the 
partial sum operator Sn from which A2m is obtained. The proof of such a sharp 
estimate is rather involved and requires detail knowledge of the zeros and weights 
of the quadrature and cubature formulas. On the other hand, a result in [22] shows 
that the norm of any projection operator from C{B'^) to U'^ is at least 0{n^~) for 
d> 2. As An in (4.2) is in fact a projection operator, its norm cannot be bounded. 
We expect that the norm is in the order of 0{n'^^'^) multiplied by a log factor. 

It should be mentioned that other polynomial based algorithms may have better 
approximation property ([15, 16]. However, the polynomial preserving property 
seems to be an important characteristic of OPED and using the partial sum allows 
also fast implementation of the algorithm. The numerical tests show that OPED 
works very well even for step functions such as Logan-Sheff head phantom [28, 29] . 

5. Singular value decomposition of the Radon transform 

Let A : H h-^ K he a. linear continuous operator, where H and K are Hilbert 
spaces. Let {fk}k>o and {gk}k>o be orthonormal systems with respect to the inner 
product (•, ■)h in H and (•, ■)k in K, respectively. The singular value decomposition 
of ^ is a representation 

oo 

(5.1) Af = ^-ik{f,fk)H9k, 

k=l 
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where are the singular values of A. Let A* be the adjoint of A. Then 

oo 

(5.2) A*g = ^lk{g,gk)Kfk- 

k=l 

Evidently Afk = 'jkQk and A*gk = 'Jkfk- Furthermore, the generalized inverse of A 
is given by 

oo 

(5.3) A+9 = J2%\f,fk)Hgk. 

The singular value decomposition of the Radon transform was developed in [5, 
12] (see also [19]). Let Z = S^''^ x [-1,1] and w{t) = Vl - 1^, and denote by 
L^{Z, w^~''') the space of Lebesgue integrable functions 

L\Z,w'-'') := {g : g{~t-t) = gi^t), \\g\\L2(z) < oo}, 

where ||ff||^2(2) = {9,9)l^{z) and the inner product is defined by 

{f,9)LHz)~Cd/2 f <Jd^ I f{^,t)g{^,t)doj{Oil-f)'^dt, 

J-i Js-i-^ 

in which 0^/2 is defined as in (2.2). Then it is known (see, for example, [19]) that 

is continuous. An orthonormal basis of L'^{Z,w^~'^) is readily available. 

Proposition 5.1. Let {Yj,™ : 1 < j < dimHJ^} denote an orthogomal basis ofUf^ 
and define 

(5.4) gU^,t) = [hW^^]-'/^(l-t')^c'J\t)Yj,r,_,k{0, 

where hlf^"^^ is defined in (2.2). Then the functions {gkj '■ < 2k < n,l < j < 
dimW^_2;.} forms an orthogonomral basis for L'^{Z,w^~'^). 

Proof. It is straightforward to verify that {,9" j} form an orthonormal system of 
L'^{Z,w^~'^). Let g G L'^{Z,w^~'^). Then w-^^^^g can be expanded in terms of 
the product orthonomal basis {[hi^^'^^]~^^'^Cn^'^(t)Yj^n^m{^) : < m < n, < 
j < dimni_.^} of L'^{Z,w'^-'^). The condition g{-^, -t) = g{^,t) shows that the 
coefficients of the expansion are zero whenever m is odd, so that we can assume 
TO = 2k and the expansion is uniquely determined. □ 

Using fj}j in (2.4) and (5-4), the singular value decomposition of the Radon 
transform at (5.1), (5.2) and (5.3) become the following: 

Theorem 5.2. Assume f is in the Schwartz space. The singular decomposition of 
TZf is 

(5.5) 7e/ = ^7n E J2 <f'fk,j)LHB^)9kj 

71=0 0<2A;<n j=0 

where Mm = dimHm, Cd/2 is defined at (2.2); and 

00 M-n -2k 

(5.6) 7^*5 = E^» E E i9,9h)LHz)fk,y 

71=0 0<2k<n j=0 
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Furthermore, 

(5.7) /(a;) = ^7-1 ^ J2 ^9 , 9h) i^Hz) fk,j ■ 

n=a a<2k<n j=0 

These equations are the reahzation of (5.1), (5.2) and (5.3) for the Radon trans- 
form. They are exactly the SVD derived in [5, 12], once the difference in notations 

is accounted for. 

Below we derive the singular value decomposition using our notation here. We 
need a proposition that goes back to [17] when d = 2. 

Proposition 5.3. Let P e V^. Then for each t e [-1, 1] and £, e 

In particular, the above formula applies to harmonic polynomials of degree n. 
Proof. Let be an orthogonal matrix whose first row is ^. Then (3.1) shows that 

nP{^,t) = {l-t^)^ [ P{{t,Vr~^y)Qi)dy. 

The integral is a polynomial of t since an odd power of Vl — t"^ is always com- 
panioned by with \a\ being odd, which has integral zero. Therefore, g{t) = 
(1 — t'^)~^~TZf{^,t) is of degree k in t. Furthermore, the integral shows that 

5(1) = vol(B''-i)P(0 = 6d-iP(0- 

If gj e for < j < n - 1, then the equation (3.2) and the fact that P e lead 
to 



(5.8) np{^,t) = ba-Ai - t-Y-^%J^p{i). 



1 r 

g{t)gj{t)il-t^)^dt= / P{x)gj{{x,O)dx = 0, 
1 Jb'I 



which shows immediately that the polynomial g{t) is an orthogonal polynomial 
with respect to (1 — t'^)^ on [—1, 1], that is, 

9{t) = {i-tY'^nfiU) = aCi/^{t). 

Setting t = 1 determines the constant a and completes the proof. Finally, (2.4) 
with A; = show that harmonic polynomials of degree n are in V^. □ 

Corollary 5.4. Let f^j be the orthonormal basis ofV^ given in (2.4). Then 

7^/fe,,(e,t)=7n5fe,,(^,i), 
where the singular values 7„ of TZf are given by 



(5.9) 7„ = bd-Wn\/{d)n. 

Proof. Using (3.5) and (5.4), the equation (5.8) shows 

(^n (1) 

where 7„ = bd-i[hn^'^^]~^/'^Hn/Cn^'^\^), which can be simplified by using (2.2), 
(3.6) and the fact that ci^^'^\l) = {d)n/n\. □ 
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Theorem 5.5. The singular decomposition ofTZf satisfies 

oo — 1 r 

(5.10) nf = c,/2{l-s')'^^\hW'A' / f{x)C^/'i{x,0)dxC^/'{t), 
where Mm = dimWJ^, C4/2 is defined at (2.2); and 

(5.11) 7^*5 = c^7nK/'^ / / 7^/(^,^)c;^/2((x,0)c;?/'(^)dc^(0^^^. 

where c = 0^/2^4-1(^4^ ■ 

Proof. To prove (5.10), we note that by (3.5), 

(5.12) fl^{x)gl^{^,t) = [hW'Y'^'H-\l - t')'^Ct'\t)fl^{x)fl^{0. 

Since the constants are independent of k and j, we can use (2.7) to write the 
summations in k and j of (5.10) in a compact form. Collecting constants and using 
(3.6), (5.9) and (2.2), we easily verify that 

^^[hi^l^)]-y-H-^fl^{x)'^±^ = 

Finally we note that bd-ib^^ = 0^/2- The proof of (5.11) is similar. □ 

It is worth to comment that the two expressions (5.10) and (5.11) arc independent 
of the choice of orthonormal bases, and the equation (2.7) implies that we can 
deduce the SVD from them using any orthonormal basis. In [5, 12], the SVD in 
terms of orthogonal basis with respect to is derived. In these more general 
cases, however, the simple analogue of the second equations of (5.5) and (5.6) do 
not hold. The reason again lies in the fact that the constant in (3.9) depends on k. 

Finally, by (5.3), the truncation of the expansion of / becomes 

m=0 0<2fc<m j=0 

Just as in the equations (5.5) and (5.6), we can use (5.12) and (2.7) to derive a 
compact formula. The formula, however, is exactly Snf- As a consequence, we see 
the truncated SVD algorithm agrees with that formula (3.11). Hence, truncated 
SVD can be effectively implemented by using the OPED algorithm. 
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